1 IMPORT

1.1 Import librairies

library(ggplot2)
library(cowplot)
library(tidyverse)
library(dplyr)
library(Seurat)

setwd("~/Projects/HumanThymusProject/")
source("~/Projects/HumanThymusProject/scripts/colors_universal.R")

1.2 Import data

# import seurat objects
seur.nkt  <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_inkt.rds")
seur.mait <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_mait.rds")
seur.gdt  <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_gdt.rds")
seur.cd4  <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_cd4.rds")
seur.cd8  <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_cd8.rds")

# sanity check seurat objects correspond to Fig 2
DimPlot(seur.nkt,  group.by="clusters_per_lineage", label=T, cols=cols_thym_nkt)

DimPlot(seur.mait, group.by="clusters_per_lineage", label=T, cols=cols_thym_mait)

DimPlot(seur.gdt,  group.by="clusters_per_lineage", label=T, cols=cols_thym_gdt)

# import GEP usage dataframe
gep_usage <- read.table("./data_github/cNMF/non_imputed_cNMF_allcells.usages.k_12.dt_0_02.consensus.txt", header=T)
colnames(gep_usage) <- paste0("GEP", c(2,5,3,1,4,12,6,7,8,10,9,11), "_usage")
gep_usage <- gep_usage[,paste0("GEP", 1:12, "_usage")]
head(gep_usage)
##                     GEP1_usage GEP2_usage GEP3_usage GEP4_usage  GEP5_usage
## CD4_1_Thymus_878524  0.0000000 0.00000000 0.01697069 0.02701780 0.003293012
## CD4_1_Thymus_679405  0.0000000 0.00000000 0.00000000 0.00000000 0.000000000
## CD4_1_Thymus_281188  0.0000000 0.01499587 0.02538062 0.01897941 0.000000000
## CD4_1_Thymus_285601  0.3893558 0.19701648 0.01374016 0.03301096 0.009603295
## CD4_1_Thymus_174537  0.0000000 0.00000000 0.01813514 0.02750585 0.001959355
## CD4_1_Thymus_128406  0.0000000 0.02029776 0.02094953 0.01324090 0.001427132
##                     GEP6_usage   GEP7_usage GEP8_usage GEP9_usage GEP10_usage
## CD4_1_Thymus_878524          0 0.0000000000 0.01670451 0.06574529  0.14537527
## CD4_1_Thymus_679405          0 0.0000000000 0.00000000 0.00000000  0.00000000
## CD4_1_Thymus_281188          0 0.0002823688 0.01075830 0.09749147  0.21138272
## CD4_1_Thymus_285601          0 0.0178349300 0.22377394 0.05961335  0.01831705
## CD4_1_Thymus_174537          0 0.0000000000 0.01610514 0.06646363  0.15172973
## CD4_1_Thymus_128406          0 0.0055636405 0.01326269 0.08876818  0.21684438
##                     GEP11_usage GEP12_usage
## CD4_1_Thymus_878524 0.655783300  0.00000000
## CD4_1_Thymus_679405 1.147190900  0.00000000
## CD4_1_Thymus_281188 0.563456800  0.00000000
## CD4_1_Thymus_285601 0.008145613  0.03549316
## CD4_1_Thymus_174537 0.649509700  0.00000000
## CD4_1_Thymus_128406 0.565640400  0.00000000

2 ANALYSIS

2.1 Add GEP1-GEP11 usage to seurat objects

# check that rownames are in same order
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.cd4@meta.data),])==rownames(seur.cd4@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.cd8@meta.data),])==rownames(seur.cd8@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.nkt@meta.data),])==rownames(seur.nkt@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.mait@meta.data),])==rownames(seur.mait@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.gdt@meta.data),])==rownames(seur.gdt@meta.data), useNA="ifany") # all TRUE

# add GEP usage to seurat objects
gep_colnames <- paste0("GEP", 1:11, "_usage")
seur.cd4@meta.data[,gep_colnames]  <- gep_usage[rownames(gep_usage) %in% rownames(seur.cd4@meta.data), gep_colnames]
seur.cd8@meta.data[,gep_colnames]  <- gep_usage[rownames(gep_usage) %in% rownames(seur.cd8@meta.data), gep_colnames]
seur.nkt@meta.data[,gep_colnames]  <- gep_usage[rownames(gep_usage) %in% rownames(seur.nkt@meta.data), gep_colnames]
seur.mait@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.mait@meta.data), gep_colnames]
seur.gdt@meta.data[,gep_colnames]  <- gep_usage[rownames(gep_usage) %in% rownames(seur.gdt@meta.data), gep_colnames]

2.2 Plot GEP usage per thymic lineage

# Plot
seur_vector <- list("iNKT"=seur.nkt, "MAIT"=seur.mait, "GDT"=seur.gdt, "CD4"=seur.cd4, "CD8"=seur.cd8)

plist_massive <- list()
for(gep in gep_colnames){
  print(gep)
  # get row of plots (for one GEP)
  plist <- list()
  for(i_seur in names(seur_vector)){
    print(i_seur)
    seur <- seur_vector[[i_seur]]
    # get legend once
    if(gep=="GEP1_usage" & i_seur=="iNKT"){
      plegend <- ggpubr::get_legend(
        SCpubr::do_FeaturePlot(seur,  features=gep, ncol=1, legend.position="right", legend.title="GEP usage", border.size=1, pt.size=2, order=T)+
          scale_color_viridis_c(limits=c(0,1.25), option="B")+
          theme(plot.background = element_rect(color = "black"))
      )
    }
    # make plot
    legendpos <- "none"
    plt_title <- ""
    if(gep=="GEP1_usage"){plt_title=i_seur}
    p <- SCpubr::do_FeaturePlot(seur,  features=gep, ncol=1, legend.position=legendpos, border.size=1, pt.size=2, order=T, plot.title=plt_title)+
      scale_color_viridis_c(limits=c(0,1.25), option="B")+
      theme(plot.background = element_rect(color = "black", linewidth = 2),
            plot.margin=unit(c(0,5.5,0,5.5), "points"))
    plist[[i_seur]] <- ggrastr::rasterise(p, layers="Point", dpi=300)
  }
  prow_gep <- plot_grid(plotlist=plist, nrow=1, scale=0.9)
  plist_massive[[gep]] <- prow_gep 
}
## [1] "GEP1_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP2_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP3_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP4_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP5_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP6_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP7_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP8_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP9_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP10_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP11_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
plot_grid(
  plot_grid(plotlist=plist_massive, ncol=1, align="h"),
  ggpubr::as_ggplot(plegend), ncol=2, rel_widths = c(5, .5), scale=c(0.95, 5)
)

3 SESSION INFO

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.2.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Paris
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] Seurat_5.0.1       SeuratObject_5.0.1 sp_2.1-2           lubridate_1.9.3   
##  [5] forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2       
##  [9] readr_2.1.5        tidyr_1.3.0        tibble_3.2.1       tidyverse_2.0.0   
## [13] cowplot_1.1.2      ggplot2_3.4.4     
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3     rstudioapi_0.15.0      jsonlite_1.8.8        
##   [4] magrittr_2.0.3         ggbeeswarm_0.7.2       spatstat.utils_3.0-4  
##   [7] farver_2.1.1           rmarkdown_2.25         vctrs_0.6.5           
##  [10] ROCR_1.0-11            Cairo_1.6-2            spatstat.explore_3.2-5
##  [13] rstatix_0.7.2          htmltools_0.5.7        broom_1.0.5           
##  [16] sass_0.4.8             sctransform_0.4.1      parallelly_1.36.0     
##  [19] KernSmooth_2.23-22     bslib_0.6.1            htmlwidgets_1.6.4     
##  [22] ica_1.0-3              plyr_1.8.9             plotly_4.10.3         
##  [25] zoo_1.8-12             cachem_1.0.8           igraph_1.6.0          
##  [28] mime_0.12              lifecycle_1.0.4        pkgconfig_2.0.3       
##  [31] Matrix_1.6-5           R6_2.5.1               fastmap_1.1.1         
##  [34] fitdistrplus_1.1-11    future_1.33.1          shiny_1.8.0           
##  [37] digest_0.6.34          colorspace_2.1-0       patchwork_1.2.0       
##  [40] tensor_1.5             RSpectra_0.16-1        irlba_2.3.5.1         
##  [43] ggpubr_0.6.0           labeling_0.4.3         progressr_0.14.0      
##  [46] fansi_1.0.6            spatstat.sparse_3.0-3  timechange_0.2.0      
##  [49] httr_1.4.7             polyclip_1.10-6        abind_1.4-5           
##  [52] compiler_4.3.1         withr_2.5.2            backports_1.4.1       
##  [55] viridis_0.6.4          carData_3.0-5          fastDummies_1.7.3     
##  [58] highr_0.10             ggsignif_0.6.4         MASS_7.3-60           
##  [61] tools_4.3.1            vipor_0.4.7            lmtest_0.9-40         
##  [64] beeswarm_0.4.0         httpuv_1.6.13          future.apply_1.11.1   
##  [67] goftest_1.2-3          glue_1.7.0             nlme_3.1-164          
##  [70] promises_1.2.1         grid_4.3.1             Rtsne_0.17            
##  [73] cluster_2.1.6          reshape2_1.4.4         generics_0.1.3        
##  [76] gtable_0.3.4           spatstat.data_3.0-3    tzdb_0.4.0            
##  [79] data.table_1.14.10     hms_1.1.3              car_3.1-2             
##  [82] utf8_1.2.4             spatstat.geom_3.2-7    RcppAnnoy_0.0.21      
##  [85] ggrepel_0.9.5          RANN_2.6.1             pillar_1.9.0          
##  [88] spam_2.10-0            RcppHNSW_0.5.0         later_1.3.2           
##  [91] splines_4.3.1          lattice_0.22-5         survival_3.5-7        
##  [94] deldir_2.0-2           tidyselect_1.2.0       miniUI_0.1.1.1        
##  [97] pbapply_1.7-2          knitr_1.45             gridExtra_2.3         
## [100] SCpubr_2.0.2           scattermore_1.2        xfun_0.41             
## [103] matrixStats_1.2.0      stringi_1.8.3          lazyeval_0.2.2        
## [106] yaml_2.3.8             evaluate_0.23          codetools_0.2-19      
## [109] cli_3.6.2              uwot_0.1.16            xtable_1.8-4          
## [112] reticulate_1.34.0      munsell_0.5.0          jquerylib_0.1.4       
## [115] Rcpp_1.0.12            globals_0.16.2         spatstat.random_3.2-2 
## [118] png_0.1-8              ggrastr_1.0.2          parallel_4.3.1        
## [121] ellipsis_0.3.2         assertthat_0.2.1       dotCall64_1.1-1       
## [124] listenv_0.9.0          viridisLite_0.4.2      scales_1.3.0          
## [127] ggridges_0.5.5         leiden_0.4.3.1         rlang_1.1.3